Many physical and engineering problems can be effectively described by non-linear symbolic models. However, existing non-linear symbolic regression (SR) methods, such as Equation Learner and Kolmogorov-Arnold Network, are constrained to a limited set of elementary continuous functions. Here, we introduce the Neural Operator-based Model Approximation and Discovery (NOMAD) method, a novel approach to symbolic model discovery that leverages Neural Operators (NOs) to approximate broad range of symbolic operations. NOMAD is structured as a computational graph that spans all potential non-linear combinations of symbolic operations, up to a specified depth of the computation graph. By learning a sparse set of coefficients, NOMAD constructs compact and interpretable symbolic expressions. We demonstrate that the proposed algorithm can discover symbolic expressions containing elementary functions with singularities, special functions, and derivatives. Furthermore, we demonstrate that NOMAD can successfully rediscover the governing partial differential equation of a two-dimensional diffusion system, showcasing its ability to handle multi-dimensional differential operators. By broadening the scope of symbolic operations available for discovery, NOMAD advances the capabilities of existing SR methods and provides a powerful and flexible tool for model discovery, capable of capturing complex dynamics in a variety of physical systems.